function [t,x,maxdist,maxheight,Vburnout] = runSimDeltaVIsp(Isp,plotThings,ispBurnoutFlag)
% function to run water bottle rocket simulation
% inputs = [thrust-data-filename, plot-decision]
% ouputs = [t,x,maxdist,maxheight,loadtime,simTime] 
% thrust-data-filename = filename corresponding to desired static thrust
% data
% plot-decision = 1: plot all results of child functions,  0: no plotting
% t = vector of simulation time 
% x = array of rocket states corresponding to simulation time state
% variables are [x-position, y-position, z-position, x-velocity,
% y-velocity,z-velocity, water-mass]
% maxdist - distance of rocket on ground impact
% maxheight - apogee height of rocket
% Vburnout, max velocity of flight, should correspond to burnout velocity


% create data structure that contains rocket and atmospheric condition
% parameters
data = loadPrelimConditions(0,plotThings); 


%calculate initial velocity from isp
ve = Isp*norm(data.gavec);
deltaV = ve*log((data.bottleMass + data.initWMass)/data.bottleMass);



v0 = data.initH*deltaV;

% initialize state variables with zero position and velocity and the
% bottle's initial water mass defined in the data structure
x0 = [0,0,0,v0(1),v0(2),v0(3),0]';

% set options for ode integrator
options = odeset('RelTol',1e-9);


%% -----------------------------------------------------------------------
% perform simulation 

% perform integration either using ode45 or forward euler (one or the other
% commented out), FE provided for simplicity in troubleshooting and
% comparison
[t,x] = ode45(@rocketODE,[0,6],x0,options,data);
% [t,x] = euler(@rocketODE,[0:.001:4],x0,options,data);

avgdt = mean(diff(t)); %calculate an average time step for the simulation


%%% --------------------- calc vburnout ---------------------------------
% if ispBurnoutFlag set to 2 or exists, use the velocity at t=0.25, else
% use the deltaV
if nargin < 3
    ispBurnoutFlag = 1;
end
    
if (ispBurnoutFlag == 2)
    % find data point corresponding to 1/4 second
    [terr,qInd] = min(abs(t-.25));
    if terr > .02
        fprintf('\n Large difference between chosen t val and 0.25 sec\n desired val for vburnout \n %f sec difference\n',terr);
    end
    Vburnout  = norm(x(qInd,4:6));
else
    Vburnout = deltaV;
end


%% -----------------------------------------------------------------------
% calculate max distance, height and other parameters, also plot results if
% desired

% find the index corresponding to ~10% through the simulation (time wise),
% used to elimate portions of the rocket path where it is near the pad so
% that the closest portion of the path to the ground can be safely used as
% the point for maximum distance (where the rocket is impacting the ground)
frontOffset = find(0.1*t(end) > t,1,'last');

% consider only the last 90% of the data to find where the rocket hits the
% ground
[dum,ind] = min(abs(x(frontOffset:end,2)));
ind = ind +frontOffset; %correct the index so that it includes the all data so it can be used easily

% calculate maximium distance based on index found above, as well as max
% height
maxdist = sqrt(x(ind,1)^2+x(ind,3)^2);
maxdistft = maxdist*3.28;
maxheight = max(x(:,2));
maxHft = maxheight*3.28;

% if plots are desired, plot the trajectory in 2d
if plotThings

    figure
    plot3(x(:,1),x(:,3),x(:,2))
    hold on
    plot3(x(ind,1),x(ind,3),x(ind,2),'*','MarkerSize',10,'LineWidth',1)
    xlabel('xpos')
    ylabel('zpos')
    zlabel('ypos')
    title('rocket trajectory')
    grid on
    axis equal
    xmax = 1.2*maxdist + 1;
    ymin = 10*min(x(:,3)) - 1;
    ymax = 10*max(x(:,3)) + 1;
    zmax = 1.2*max(x(:,2));
    axis([0,xmax,ymin,ymax,0,zmax])
    view([0,0])
    
    speed = sqrt(x(:,4).^2 + x(:,5).^2+x(:,6).^2);
    figure, plot(t,speed)
    xlabel('t'), ylabel('speed')
    
    figure, plot(t,x(:,4),t,x(:,5),t,x(:,6)),xlabel('t'),ylabel('velocity'),legend('x','y','z')
    grid on
    
    
end

% if plots are desired also print the cpu times as well as max height and
% distance
if plotThings
    fprintf('The Average time step = %2.1e s\n',avgdt)
    fprintf('The max distance is %3.1f ft or %3.1f m\n',maxdistft,maxdist)
    fprintf('The max height is %3.1f ft or %3.1f m\n',maxHft,maxheight)
    fprintf('The burnout velocity is %3.1f m/s \n',Vburnout);

end